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Abstract 

This paper deals with the numerical modeling of transient mechanical waves 
in linear viscoelastic solids. Dissipation mechanisms are described using the 
generalized Zener model. No time convolutions are required thanks to the 
introduction of memory variables that satisfy local-in-time differential equa- 
tions. By appropriately choosing the relaxation parameters, it is possible 
to accurately describe a large range of materials, such as solids with con- 
stant quality factors. The evolution equations satisfied by the velocity, the 
stress, and the memory variables are written in the form of a first-order sys- 
tem of PDEs with a source term. This system is solved by splitting it into 
two parts: the propagative part is discretized explicitly, using a fourth-order 
ADER scheme on a Cartesian grid, and the diffusive part is then solved ex- 
actly. Jump conditions along the interfaces are discretized by applying an 
immersed interface method. Numerical experiments of wave propagation in 
viscoelastic and fluid media show the efficiency of this numerical modeling for 
dealing with challenging problems, such as multiple scattering configurations. 
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1. Introduction 



Wave motion in real media differs in many aspects from motion in an ide- 
alized elastic medium. The dispersion and attenuation induced, for instance, 
by grain-to-grain friction can greatly affect the amplitude of the waves and 
their arrival times. In the case of small perturbations, linear models of vis- 
coelasticity provide reasonably accurate means of describing these effects. 
Viscoelastic constitutive laws give the stress in terms of the past strain rate 
history. 

The aim of this paper is to simulate the propagation and the diffraction 
of transient viscoelastic waves. We propose a new approach in three steps: 

(i) The generalized Zener model is used fsj]. The convolution products are 

then replaced by a set of local-in-time differential equations coupled 
with the evolution equations of velocity and stress. Moreover, usual 
attenuation laws can be approximated closely. 

(ii) The evolution equations are splitted into two parts: a propagative part, 

which is solved using a fourth-order finite-difference ADER scheme on 



a Cartesian grid 33|; and a diffusive part, which is solved analytically. 



Doing so ensures an optimal condition of stability. 

(iii) The jump conditions along the interfaces are discretized by an immersed 
interface method, which introduces a subcell resolution of the geometry 
and maintains the convergence rate of the scheme despite the non- 
smoothness of the solution. See ^3] for an overview of these methods. 

The generalized Zener model (i) has been addressed by various means, such 
as finite difference methods 31, 35|, spectral methods [3], spectral-element 



methods [loj, finite element methods [IgI, 0], to cite only a few. The steps 
(ii) and (iii) combine the computational efficiency of Cartesian grid meth- 
ods and an accurate description of the interfaces, as stated in the case of 



non-dissipative media |25l . |26| and applied to computationally challenging 



configurations 10 



The article is organized as follows. In section O the generalized Zener 
model is presented; the method used to determine its parameters to sim- 
ulate a given quality factor is also described. In section |3l the evolution 
equations are written in the form of a first-order hyperbolic system with a 
source term; the jump conditions along the interfaces are also stated. The 
numerical methods used are introduced in section Hj including the numerical 
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scheme and the sphtting for the integration of the evolution equations, and 
the immersed interface method for the discretization of the jump conditions. 
Numerical experiments are presented in section [5] in the case of a viscoelastic 
/ fluid interface. Comparisons with analytical solutions are proposed. A nu- 
merical experiment involving multiple scattering in a random medium also 
confirms the efficiency of the approach. Lastly, the perspectives are discussed 
in section [HI 

2. Physical modeling 

2.1. Constitutive law 

In a viscoelastic solid undergoing small perturbations, the stress depends 
linearly on the history of the past strain rates. In ID, one writes 



where a is the stress, e = ^ is the strain, u is the displacement, ip(t) is the 
relaxation function, and * denotes the time convolution. 

Various models of viscoelasticity can be found in the literature jsf. The 
Maxwell model predicts a vanishing asymptotic residual stress. It therefore 
appears more appropriate for representing viscoelastic fluids. The Kelvin- 
Voigt model is computationally advantageous joj, but it predicts an un- 
bounded phase velocity as frequency increases. Here, we choose the gen- 
eralized Zener model, which accurately mimics the mechanical behavior of 
classical viscoelastic media during relaxation experiments: 



where H refers to the Heaviside distribution, Nr is the number of relaxation 
mechanisms, di are relaxation frequencies, and the coefficients Ke are strictly 
positive. The instantaneous unrelaxed state is denoted by ir^, and at the end 
of the process, the relaxation function has returned completely to the positive 
relaxed modulus vr^, where < nr < Hu- The phase velocity increases with 
the frequency, from cq = \/viv7p ^^^^ frequency to Cqo = ^/tTu/p at infinite 
frequency, where p is the density 




(1) 
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2.2. Determination of the parameters 

Let J-" be the Fourier transform of a function g{t) 



+00 

iu}t 



F{g)= / g{t)e-'^Ut, (3) 



00 



where u is the angular frequency. From ([2]), the modulus of viscoelasticity 
M(u) = is: 

MM^..(i+..|:^). (4) 

and the ratio between the imaginary and real parts of M is: 

1^ i^e 



The quality factor Q characterizes the attenuation of the viscoelastic waves. 

To determine the 2 Nr coefficients ne and 6'^ in ([2]), we choose to minimize 
the distance between Q~^{u) and a given ^ band of angular fre- 

quencies [uq, cui]. Here we have implemented a classical linear least-squares 
minimization procedure in the L2 norm 13, Q, 0|- Relaxation frequencies 
are distributed linearly on a lo gari thmic scale of Nr points, ranging from 
/o = a;o/(2vr) to /i = / (2 vr) 
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The angular frequencies loq and cui obviously depend on the spectra of the 
source. The coefficients ki are then obtained by solving the over-determined 
linear system deduced from (jS]) 

{Oe-UlkQrefi^k)) n-l f ~ \ U ^ OAT -\ 
/)2 I ~2 f^l = Qref{^k), /c = 1 , . . . , 2 A^^, - 1 , (7) 
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where Uk are distributed linearly on a logarithmic scale oi 2 — 1 points 



u:k = ujo[-^] , fc = 1,..., 2iV,-l. (8) 

In our numerical experiments, calculations based on ([7]) have never yielded 
non-physical negative values of k^, even with highly attenuating media (typ- 
ically Qref = 5). If necessary, more sophisticated methods can be applied. 
For instance, a nonlinear least-squares constraint optimization was used in 
jil to ensure that the coefficients ki were positive. An alternative method of 
optimization in the norm Lqo was presented in j2|, but this latter method is 
restricted to materials with constant Qref- 

To determine the number of relaxation mechanisms A^r, one can compare 
Qref and the quality factor Q deduced from ([5]) after optimization. However, 
the modeling error in the time domain is not easily deduced. A second idea 
consists in comparing the transient 1-D analytical solutions associated with 
Qref and Q. These solutions are calculated using classical Fourier techniques, 
which are not described here. 

2.3. Numerical examples 

The determination of the parameters in ([2]) is illustrated in figure [H The 
set up is the same here as in section |5l the source is a smoothly truncated 
sinusoid of central frequency fc = 50 Hz; the optimization of the quality 
factor is done between /o = /c / 10 and /i = 10 fc, the physical parameters 
are p = 1200 kg/m^, cq = 2800 m/s, Qref = 20. Considering a constant 
Qref is usual in geosciences, where the real media show a quasi-constant 
quality factor within very wide frequency ranges p^]. In addition, the exact 
solution associated with a constant Q^e/ is particularly simple to obtain. 



involving the Kjartansson formula |20|. Note that non-constant Qref can 



also be considered in numerical modeling without any restrictions. 

Figure [1] compares the quality factor (a) and the time-domain exact so- 
lutions (b), when A^^ = 1 or A^^ = 3. Increasing A^ clearly decreases the 
error introduced by describing the constant quality factor medium in terms 
of a finite number of relaxation mechanisms. At the same time, increasing 
Nr greatly increases the computational cost, especially the memory require- 
ments. In practice, A^^ > 3 is rarely implemented in the literature, and 
A^r- = 1 is widely used, especially in the 3-D context (sif . 

Dispersion and relaxation curves shown in figure [1] (c) and (d) are com- 
puted taking A^^ = 3 relaxation mechanisms. One observes the predicted 
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Figure 1: Determination of viscoelastic parameters, with Qref = 20, fc = 50 Hz, and 
A'r = 1 or iVr = 3 relaxation mechanisms. Quahty factors (a); the sohd horizontal line 
gives the exact value 1 / Qref- Time-domain 1-D analytical solutions obtained with the 
constant-Q model (Kjartansson's model) and the generalized Zener model (b). Phase 
velocity (c) and relaxation function (d) obtained with Nr = 3. Physical parameters 
are those used in section [5l In (a) and (c), the dotted vertical lines give the relaxation 
frequencies when Nr — 3. 
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behavior: strict increase of the phase velocity from cq to c^o (c), strict de- 
crease of the relaxation function from vr^ to tt^ (d). 



3. Initial boundary-value problem 

3.1. Constitutive law in two- dimensions 

The viscoelastic law ([1]) is generalized so that it applies to all space di- 
mensions. In the 2-D case, the constitutive law governing a linear isotropic 
viscoelastic medium is ISjl 



cr. 



mt) -2Mt)) * ^'5.. + 2Mt) * (9) 



where aij and the components of the stress and strain tensors, and Sij 

is the Kronecker symbol. With the generalized Zener model, the relaxation 
functions ipn and are given by 



Nr 



Mt)=f^r(l + f^Kte~'A Hit) 



^=1 / (10) 



where tt,. = pc^^ and pr = pc^^ are relaxed moduli under compressional and 
shear loads. The phase velocities of the compressional (P) and shear (S) 
waves at zero frequency are denoted Cpo and Csq. The unrelaxed moduli are 
written 




Nr 



7r„ = vrJ 1 + = pcj^, /i„ = /iJ 1 + = pc^^, (11) 



£=1 



where Cp^ and Cg^ are the phase velocities of P and S waves at infinite 
frequency. The parameters 9i, k,^ and in (ITU]) are determined as in section 
12. II from the quality factors Q-^^j and Q^^j of P and S waves. Usually, Q^^y < 
Q^ef- ^ waves are more attenuated than the P waves. The relaxation 
frequencies Og are the same with both P and S waves, since they depend only 
on the frequency band of interest ([6]). In addition, describing P and S waves 
with identical relaxation times, as well as identical numbers of relaxation 



mechanisms, greatly reduces the memory requirements [SlJ, |35 
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3.2. Evolution equations 

To obtain the evolution equations satisfied by (Tjj, the constitutive law 
([9]) is differentiated in terms of t, taking (|T0|) . If z = j, we obtain 



7r„ - 2 /i„ - — + 2 /i„ -— + > ^ij-^, 12 



t O Xk O Xj 



where the ^ije are called memory variables 
6i^ = -^^ (7r,4-2/i,«:|)e-^^*i/(t)*-^-2/i,0,/€|e-^^*ff(t)*-^. (13) 

OXk OXj 

These memory variables satisfy the differential equations 

- (TTr^ h2/i^K^- , i = l,...,Nr. 



dt ^ V dXk 'dx, 

In the same way, if i 7^ j, we obtain 



(14) 



dan f dvi dvj\ v-^ ^ 



dt \dxj dxi. 



with the memory variables 

6,, = -fir 9, 4 e-'^ 'H{t)*(^ + P-), (16) 

\OXj OXiJ 

that satisfy the differential equations 

- -e, fe.,. + ^^r4(^ + p-\], i = 1, AT,. (17) 



dt ^ V V-^a;, 



The convolutions in (fT3!) and (fT6|) induced by the convolution in ([9]) are 
no longer involved in f[T^ and (1171) : adding a set of memory variables that 
satisfy local-in-time differential equations avoids to store the past values of 
the solution. In 2-D contexts, combining (|T2l) . (|T^ . f lTSj) and f|T71) with 
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Newton's law yields a system of 5 + 3 A^^ partial differential equations 
( dvi I {dan ^ dai2 _ ^ 



dt p \ dx dy 
dv2 1 f dai2 ^ da22 , _ ^ 



dt p \ dx dy 
dan dvi ^ .dv2 

dt dx dy ^-^ 

^ 1=1 

dai2 ( dvi dv2\ ^ 

/ dvi dv2 sr^ . 

dt dx dy ^-^ 

— — + Oe TX,. K^— h [iXr K'^-^Pr Kg) — — = -9e , i=l, Nr 

dt \ dx dy J 

9^121 . a s f 9vi dv2\ 
dt \dy dx J 

d^22i , n (( P o s^dVl pdV2\ 

—— + 6'^ {iTr k'^ - 2 Pr Kf,) — h VT^ ^ ^— = -0e^22e, I,---, ^'r- 

dt \ dy J 

(18) 

Setting 

U= {Vi, V2, ail, Cri2, Cr22, ^lU,---, ^llNr, Cl21,---, il2Nr, C22I, ••• C22Afr) ) 

(19) 

one can write (fTS!) in the form of a first-order linear system with a source 
term 

— U + A — U + B — U = -SU, (20) 
dt dx dy 

where A, B and S are (5 + 3 A^^) x (5 + 3 A^^) matrices. The eigenvalues of A 
and B are real: ±Cpoo, ±Csoo, and with multiplicity 3 A^r + 1- As deduced 
from (E]), the spectral radius of S is 

^(S)=^^. = ^ = /i. (21) 

Z TX 

For further use, we introduce the restriction of U to the velocity and stress 
components and without any memory variables: 

U = {vi, V2, ail, CF12, a22f ■ (22) 
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An overline is also used to denote the restricted 5x5 matrices A and B 
involving only the velocity and stress components. 

Even in non-viscoelastic subdomains, the evolution equations are written 
in the same way as fl20|) . For instance, Qi is assumed to be a fluid medium 
in section 13.31 and in the numerical experiments. In this case, U = U = 
(fi, V2, p)'^, where p is the acoustic pressure, A and B are 3x3 matrices, 
and S = 0. Lastly, subscripts will be used to denote the medium under 
investigation: as an example, Aq is the matrix A in f2o. 

3.3. Interface conditions 




n 



Figure 2: Interface T between two media fig et ili. 



The physical parameters defined in section 13.21 can vary discontinuously 
across interfaces. In what follows, we will focus on two domains Qo and 
fli, which are separated by a stationary interface F described by a paramet- 
ric equation (x(r), ?/(t)) (figure [2]). The domain Qq contains a viscoelastic 
medium described by the generalized Zener model: the constitutive law and 
the evolution equations in Qq are described in sections 12.11 and 13. 2[ The 
domain Qi can contain vacuum, a perfect fiuid, an elastic solid or any other 
viscoelastic medium: all these combinations have been implemented numer- 
ically and tested. In the rest of the study, we will focus on the case where 
Qi contains a fiuid. In this case, the interface conditions are 

[v.n] = 0, (a.n).n = -p.n^ ((T.n).t = 0, (23) 

where [.] refers to the jump from Qq to Qi, and the unit tangential vector t 
and the unit normal vector n are 

t = = (x, y') , n = = (y, —x'^ . (24) 
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Derivatives x = ^ and y' = ^ are assumed to be continuous everywhere 
along r, and to be differentiable as many times as required. 

The interface conditions fl23|) only involve velocity and stress components 
U (|22|) . This is also the case for other configurations, such as viscoelas- 
tic / vacuum or viscoelastic / viscoelastic interfaces. In section 14.21 and 



Appendix B, we will write the interface conditions satisfied by the spatial 
derivatives of U up to the k-th order, and hence the following notation is 
introduced: 

U,^= hm fu^..., , ^" , U^..., T^U^y, (25) 

where a = 0, /c, /3 = 0, a, and £ = 0, 1 denotes the number of the 
medium Vti. The interface conditions are then written 

c;u; = c°u°, 

(26) 

L?U? = 0, L°U° = 0, 

where are the matrices of the jump conditions, and L° are the matrices 
of the boundary conditions. With this formalism, the viscoelastic / fluid 
interface conditions ( 123|) and the vectors (12^ yield 

CS(r) = 



C?(r) 




(27) 



L°(r) = (0 x'y y' - x" -x y ) , L;(r) = ( O). 

To conclude on that topic, it is noticed that no interface conditions are im- 
posed on memory variables. In some cases, however, authors exhibit bound- 
ary conditions satisfied by ^iji. see, for instance, equation (4) of j32| in the 
case of viscoelastic / vacuum interface. Such conditions are not required 
to get a well-posed problem. Moreover, they are deduced from the original 
boundary conditions satisfied by a, the constitutive law ([9]), the positivity of 
relaxation functions, and the definition of memory variables f|T3|) and (|T6|) . 
These additional boundary conditions are not useful in the immersed inter- 
face method (section | 
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4. Numerical modeling 

4.I. Numerical scheme 

Generalities. Let us take a uniform grid, with the spatial mesh size 
Ax = Ay and the time step At. An approximation U"^ of \J{xi = 
iAx, Uj = j Ay, tn = nAt) is sought. The numerical methods recalled 
in section [T] usually consist in simultaneously discretizing the propagating 
part and the source term in fl20l) . This approach has two drawbacks. First, 



building unsplit methods for ffTSj) is a difficult task [2l|, whereas large classes 
of methods already exist for hyperbolic systems without the source term S U. 
Secondly, a Von-Neumann stability analysis typically yields 

f 1 Ax 2 \ 
At<min , , (28) 



c. 



poo 



where 7 depends on the scheme. Based on fl^ and (1251) . the spectral radius 
of S induces a more restrictive bound than the classical CFL condition if 
/i > 2cp^ / (7 Ax), where /i is the maximum frequency considered during 
the determination of the parameters (section 12. 2p . The efficiency of the 
scheme is therefore penalized if large values of /i are taken. 

Splitting. Here we choose another strategy based on solving alternatively 

u L u X (J y 

g (29) 

-U = -SU. w 

The discrete operators used in stages (a) and (6) are denoted by and H^, 
respectively. The algorithm of A/'-th order splitting is written 

. ug™-^) = H,(c^At)u!7-^\ m = l,...,^^ 

( 30 ) 

Ug"^) = Il,{d^At)V^^p-'\ 

Tjn+l ^ TT(2Ar) 

or equivalently 

^IV = f n (dM-m+i At) oUa (cAr-™+i A t) \ Vl^. (31) 
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The coefficients and in fl3Up and f l3ip are given in Appendix A The 
cases M = 2 and = 4 are called Strang splitting and Ruth splitting, re- 
spectively. 



Solvers. To solve the propagative stage fl29|) -(a). many standard solvers 
for hyperbolic systems can be used as a discrete operator H^. Here we 



choose a fourth-order ADER scheme [33|. This is an explicit two time step 



spatially-centered flux- conserving scheme, with a centered stencil of 25 nodes. 
On Cartesian grids, this scheme amounts to a fourth-order Lax-Wendroff 
scheme. It is dispersive of order 4 and dissipative of order 6, with a stability 



limit 7 = max(cp^ At/Ax) = l 24 



The diffusive stage f l29p -(&) is solved exactly. Based on the notations of 
f pOj) . we obtain 



„,(2m) ^ (2m-l) 

Nr 

A2m) ^ (2m- 

Pi ""pi ' ^0 

=1 ^ 
m 

^pql ~ ^ '^pql 



<r = <r"'^ + E } (1 - ^"'^'"^^0 (32) 

f(2m) _ At f(2m-l) 



where (p, g) = {(1, 1), (1, 2), (2, 2)}, £ = 1, A^^ and m = 1, ■ ■ ■ , AT. The 
computation of e~^*'^™^* is time-consuming. To make the computational 
time of the diffusive stage (6) negligible, the exponentials in (132|) are com- 
puted and stored once in each viscoelastic subdomain and at each time sub- 
step. 

Choice of M. The splitting ( l29l) with operators Hq and H;, is A/'-th order 
accurate. If A/" = 1 or A/" = 2, the optimal CFL condition of ADER scheme 
is recovered:7 = 1. If A" = 3 or A/" = 4, then the stability limit is improved: 
numerical experiments indicate 7 ~ 1.54 and 7 ~ 1.60, respectively. To 
choose A/", the following remarks are done: 

1. for almost the same computational cost. A/" = 2 is twice accurate than 
A^ = 1, and hence Strang sphtting is preferred to first-order splitting; 

2. A/" = 4 maintains the fourth-order accuracy of ADER scheme. However, 
four spatial integrations are required, compared with a single one when 
U = 2- 

3. in counterpart, coarser spatial and temporal grids can be used with 
Ruth splitting to ensure the same accuracy than with Strang splitting: 
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the CPU time and the number of time steps are reduced accordingly. 
Moreover, the CFL hmit of stabihty is improved with Ruth sphtting. 

Further discussions on this topic and convergence measurements will be pro- 
posed in section | 



4-2. Immersed interface method 

To solve the propagative part (!29|) -(a) accurately, the solution must be 
sufficiently smooth on the whole stencil. At the irregular points, where the 
stencil crosses the interface, the smoothness requirement is not satisfied, and 
the discrete operator Hq (15U]1 needs to be modified to maintain the accuracy. 



For that purpose, an immersed interface method is implemented [30|, |25|, |26 



The basic principle of this method is as follows. Let us take an irregular 
point (xj, i/j) G Qq. The numerical computation of Uj-^™" in f pOj) requires 

to use a value at (x/, yj) G (figure[3]). Instead of using U^^J'"^'', a modified 
value J is injected into the discrete operator H^. It amounts to a k-th 
order extension of the solution from Qq into Qi, where k is an integer to be 
defined. 

In a few words, IJ} j is build as follows. Let P be the orthogonal projec- 
tion of (x/, uj) on F, and consider the disc V centered on P with a radius q 
(figure [3]). Based on the interface conditions (126|) at P and on the numerical 

values u/7 ^ at grid nodes inside P, a matrix M. is defined so that 

\JIj = M (V^^""'^^^ ^. (33) 

(2m— 2) 

It is noticed that restrictions (1221) to velocity and stress components Uj- j 



are considered in f l33|) : the system f[T8|) implies that no spatial derivatives 
are applied on the memory variables, and hence the numerical integration of 
(EH])- (a) at ( Xj, Uj) does not involve values of the memory variables at other 
nodes. 

During the propagative part fl2^ -(a). the viscoelastic medium behaves 
like an elastic medium. The derivation of the matrix M. in fl33l) is therefore 



in line with the algorithm developed for non-dissipative media [25|, with 



appropriate physical parameters. Details are given in Appendix B The only 



slight modification compared with the elastic case concerns the Beltrami- 



Michell equations: see step 2 of Appendix B 



Some comments are done about the immersed interface method: 
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Figure 3: M(x/, yj) G ili is a grid node where a modified value U| j is computed; P 
is the orthogonal projection of M onto the interface F. The grid nodes used to compute 
Uj J are inside the circle with radius q and centered on P; they are denoted by +. 



1. A similar algorithm is applied at each irregular point along F and at 
each propagative part of the splitting algorithm ( l30l) (m = 1, ■ ■ ■ , Af). 
Since the jump conditions do not vary with time, the evaluation of 
the matrices in ( l33l) is done during a preprocessing step. Only small 
matrix-vector products are therefore required at each splitting step. 
After optimization of the computer codes, this additional cost is made 
negligible, lower than 1% of the time-marching. 

2. The matrix Ai in f l33|) depends on the subcell position of P inside 
the mesh and on the jump conditions at P, involving the local geom- 
etry and the curvature of F at P. Consequently, all these insights are 
incorporated in the modified value (133|) . and hence in the scheme. 

3. The simulations indicate that the number of grid nodes inside the disc 
V has a crucial influence on the stability of the immersed interface 
method. Here we use a constant radius q. Taking k = 2, numerical 
experiments have shown that q = 3.2 Ax is a good candidate, while 
q = 4.5 A X is used when k = 3. 

4. The order k plays an important role on the accuracy of the coupling 
between the immersed interface method and a r-th order scheme. If A; > 
r, then a r-th order local truncation error is obtained at the irregular 
points. However, k = r — 1 suffices to keep the global error to the r-th 
order [f3], and hence A; = 3 is used for the ADER 4 scheme. 
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5. Numerical experiments 

5. 1 . Configuration 

Here we focus on viscoelastic / fluid configurations. Tlie pliysical param- 
eters in the viscoelastic medium VLq are: 



p = 1200 kg/m^ Cp, = 2800 m/s, = 1400 m/s, Q^,. = 20, Q^,. = 15, 



with iVj. = 3 relaxation mechanisms, and in the fluid medium fii they are: 



The time evolution of the source is given by a combination of truncated 
sinusoids 



where (3^ = 2"^'^, Uc = 27t f^, the coefficients am are: ai = 1, a2 = —21/32, 
03 = 63/768, 04 = —1/512. This source is and has a central frequency 
fc = 40 Hz. Once propagation has occured across a viscoelastic medium, 
waves emitted by the source are deformed compared with (!36l) . After op- 
timizing Kg between /o = /c/10 and /i = 10 /c (E]), then ( ITTl) yields the 
high-frequency limits Cp^ = 3161 m/s and c^^ = 1645 m/s. 

The computations are performed on x A^^^ grid nodes. The discretiza- 
tion mesh isAa; = Ay = lm, and Strang splitting is used. The time step 
follows from Cp^ At/ Ax = 0.85. A third-order immersed interface method is 
used: k = 3 and q = 3.2 Ax in section [^^21 Time-domain exact solutions are 
not available in dissipative media; they are therefore computed by Fourier 
synthesis on Nf modes, with a frequency step A /. 

On the plates, —p and an are shown in fluid and viscoelastic media, 
respectively. P waves and S waves are displayed with a green-red palette and 
a yellow-magenta palette, respectively. The distinction between these waves 
is based on numerical estimates of div v and curl v. The position of a slice 
is denoted on the plates by horizontal segments. 



(34) 



p = 1000 kg/m^ c = 1500 m/s. 



(35) 




(36) 



otherwise 
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5.2. Test 1: ID medium 

As a first experiment, wave propagation is simulated in a one- dimensional 
homogeneous viscoelastic medium [0, 400] m. The physical parameters are 
those of P- waves in The initial data is the right-going part of the field 
emitted by a source point at x = and propagated during 0.05 s (figure HJ-a). 
This field is computed by Fourier synthesis. 
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Figure 4: test 1, homogeneous viscoelastic medium. Stress a at the initial instant (a) and 
after 200 time steps (b). 

Figure HJ-b shows the stress obtained by taking A^^^ = 400 grid nodes and 
Nt = 200 time steps. One clearly observes the attenuation and the dispersion 
induced by the viscoelastic constitutive law. Excellent agreement is observed 
between numerical and exact values. 

To study more quantitatively the accuracy of the numerical scheme, con- 
vergence measurements are performed by considering two different splittings: 
Strang splitting (A/" = 2) and Ruth splitting (A/" = 4); see section WTl for 
details. Special care is taken to ensure that the initial data and the exact so- 
lution are valuable reference solutions: Nf = 65536 Fourier modes are used, 
with a frequency step A / = 0.01 Hz. Errors in norm I2 and convergence 
rates are reported in table [H and then are displayed in figure |5l For coarse 
grids {Nx < 400), theoretical orders are not yet reached. Moreover, Strang 
splitting is more accurate than Ruth splitting. For A^; > 400, second-order 
and fourth-order rates are very closely reached, and Ruth splitting becomes 
competitive. 
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Table 1: test 1, honiogeneous viscoelastic medium. Convergence measurements. 
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Figure 5: test 1, homogeneous viscoelastic medium. Error versus the number of grid nodes. 
Slopes + 2 and +4 are obtained with Strang splitting and Ruth splitting, respectively. 

The case of a ID heterogeneous fluid / viscoelastic medium is illustrated 
in figure O The interface is located at a; = 200 m. The initial right-going 
wave is put in the fluid medium (figure [6]-a) . Wave propagation is simulated 
during 200 time steps. At the final instant, the incident wave has interacted 
with the interface, and reflected and transmitted waves have been generated. 
Numerical and exact solutions are compared successfully (figure [6]-b) . 

Convergence measurements have been performed in this heterogeneous 
case, but the results are not so sharp than in table [1] and figure [51 When the 
dissipation is small {Q > 100), the results are similar to those obtained in 
the homogeneous case: second-order and fourth-order rates are obtained by 
Strang splitting and Ruth splitting, respectively. But when the dissipation 
becomes important, convergence rates between 1 and 2 are obtained, what- 
ever the splitting. The problem probably follows from the coupling between 
the splitting and the immersed interface method, during the time derivatives 



18 



of interface conditions (step 1/4 of Appendix B). These derivatives are based 
on the conservation law (12^ -(a). corresponding to a relaxed elastic medium, 
instead of the original system f lTSj) . Further analysis on that topic is required, 
in order to propose A/'-th order algorithms in heterogeneous cases. 
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Figure 6: test 1, fluid / viscoelastic medium. Stress a at the initial instant (a) and after 
200 time steps (b). 



5.3. Test 2: plane wave on a plane interface 

The second test is conducted on a 2D plane interface between the fluid 
and the viscoelastic medium. The angle between the straight line and the 
horizontal axis is equal to 70 degrees. A homogeneous acoustic plane wave 
(IP), having a wave vector inclined at an angle of 10 degrees, propagates 
in the fluid and interacts with the interface. The viscoelastic transmitted 
compressional (TP) and shear (TS) waves are inhomogeneous waves, whose 



wave vector forms a non-null angle with the direction of attenuation. See |23l . 
11 . [lil . [il, [g], i^l for further details on this topic. Figure [7] shows the incident 
field (a-c) and after 700 integration time steps (b-d), which corresponds to 
roughly 6 propagation wavelengths. Excellent agreement is observed between 
the numerical and the exact values. 

5.4. Test 3: plane wave on a circular interface 

The immersed interface method depends on the curvature of the interface 



and its successive derivatives 25, 26 . To test the method with a non-null 



curvature, we now examine a circular interface with a radius of 60 m. The 
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Figure 7: test 2. Plane interface between a fluid (on the left) and a viscoelastic medium 
(on the right). Initial instant (a-c) and after 700 time steps (b-d). IP: incident homo- 
geneous acoustic wave; RP: reflected inhomogeneous acoustic wave; TP, TS: transmitted 
inhomogeneous compressional and shear viscoelastic waves. 
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Figure 8: test 3. Circular interface between a fluid (outside) and a viscoelastic medium 
(inside). Initial field (a-b), after 220 time steps (c-d) and 440 time steps (e-f). 
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Figure 9: test 4. Initial (a) and scattered field after 2200 time steps (b). 

fluid and the viscoelastic medium are outside and inside the circle, respec- 
tively. The incident field is a plane acoustic wave with a horizontal wave 
vector. Figure IS shows the field at the initial instant (a-b), after 220 time 
steps (c-d) and after 440 time steps (e-f). Classical conversions and scat- 
tering phenomena are observed. Excellent agreement is observed between 
the numerical and the exact values. The latter are computed using Fourier 
techniques and decomposing the plane waves on the basis of Bessel functions. 

5.5. Test 4^ plane wave and multiple scattering in random medium 

In the previous tests, the validity of the numerical scheme and the im- 
mersed interface method was confirmed in the case of academic configura- 
tions. We now take a complex medium composed of 60 viscoelastic cylin- 
ders randomly embedded in water. The computations are performed on 
2000 X 3000 grid nodes. Figure [9] shows the initial field (a) and the scattered 
fields after 2200 time steps (b), when the incident wave has propagated over 
a distance corresponding to 12 wavelengths. 
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6. Perspectives 



The propagation of mechanical waves in dissipative solids was addressed 
numerically in the time domain. To avoid dealing with convolution products, 
memory variables were introduced. Evolution equations were splitted into 
two parts: the propagative part was solved numerically using a fourth-order 
scheme for hyperbolic systems; and the diffusive part was solved exactly. 
The jump conditions were discretized by means of an immersed interface 
method, which introduced a subcell resolution on a Cartesian grid. In nu- 
merical experiments, focus was put on the fluid / viscoelastic interface, but 
the algorithms have been implemented and tested in many other cases, such 
as viscoelastic / vacuum and viscoelastic / viscoelastic interfaces. 

The numerical methods presented here make it possible to simulate phys- 
ically relevant numerical experiments, for instance multiple scattering in ran- 
dom media as performed in section 15.51 By applying signal processing tools 
on the simulated data, it is possible to determine the properties of the effec- 
tive medium which is equivalent to the disordered medium investigated 10 . 
This numerical approach can be used advantageously instead of the methods 
usually adopted by physicists so far: real experiments are expensive, and 
analytical methods can be used only with very small concentrations of scat- 
terers. The latter limitation is particularly penalizing in the case of concrete, 
where the concentration of aggregates lies typically around 40%. 
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Appendix A. Coefficients of splitting 

The coefficients Cm and dm involved in (!30|) and ( 13T|) are detailed up to 
Af. They satisfy 

^ ^ Cm 1; ^ ^ dm !• 

m,=l m=l 

For = 1, one has the usual coefficients 

ci = 1, di = 1. (A.l) 
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For Af = 2, the classical second-order Strang's splitting is recovered 
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,^0, .,^1/2. 

C2 = l, C?2 = l/2. 



Since ci = in ( ]A.2p . the first spatial integration in ( 130|) is not performed. 



For = 3, the set of coefficients is 15 



ci = 7/24, di = 2/3, 

C2 = 3/4, d2 = -2/3, (A.3) 
C3 = -l/24, 4 = 1. 

Lastly, setting x =^^^ + 2"^^^ -l)/6^ 0.1756, the coefficients of fourth- 
order splitting are [15 



ci = x + l/2, rfi = 0, 
C2 = -X, = 2x+ 1, 

C3 = -X, «3 = -4x-l, 

C4 = x + l/2, 4 = 2x + l. 

Appendix B. Four steps to build Uj j ([33]) 

Step 1: high-order interface conditions. First, we seek the interface 
conditions satisfied by the spatial derivatives of the velocity and stress com- 
ponents at P. For this purpose, the zero-th order interface conditions f l26|) 
are differentiated in terms of t. The time derivatives are replaced by spatial 
derivatives, using the propagative part f l2I?]) -(a). Equations f l2B]) are also dif- 
ferentiated in terms of r, using the chain-rule. For instance, the boundary 
condition Lq Uj] = results in 

^(LSUS)^-LSA.AUS-LSB„|-US^0. 

A ,LS US) ^ (A LS) US + LS (.'A US + y§- US) ^0. 

From (IB . 1 p . a matrix Lq is built such that LJUq = 0. This matrix depends 
on r and on the physical parameters on Qq. Applying a similar procedure to 
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the three equations in fl26l) gives a set of first-order interface conditions. By 
iterating this process k times, we obtain the k-th order interface conditions 

C^Ut = CSuS, L^^U^^ = 0, £ = 0,1. (B.2) 

When k > 2, building the matrices and is a tedious task, which can 
be greatly simplified by using computer algebra tools. Note lastly that 
and involve the spatial derivatives and (a = l,..., fc + l), which 
provides insights about the local geometry of T at P. 



Step 2: high-order Beltrami equations. In ( I29l) -(a). the viscoelastic 
medium behaves like an elastic medium with Lame coefficients A = Tr^ — 2 
and fi = flu, where compatibility conditions are satisfied between some spatial 



derivatives of the strain components [28|. When expressed in terms of a, these 



conditions lead to the Beltrami equation 

Q2 Q2 Q2 



dxdy 
where 



(Ti2 = 02 



a2 



(Til + «1 



■ 0-22 + «1 



0"ll + ^2 



dy' 



0-22, 



(B.3) 



Poo 



TT,, 



■ 2/in 



4 {cl - cl ) 



(B.4) 



The equation fIB.SP is satisfied anywhere in fig. Under suitable smooth- 
ness requirements, it can be differentiated as many times as necessary, with 
respect to x and y. Since the equations thus obtained are also valid along F, 
they can be used to obtain a minimum number of independent components 

\5\ = Q\V\, £ = 0,1. (B.5) 

The algorithm for building the matrices presented in 26| can be easily 
adapted to flB.3P - flB.4p . If fii is not a viscoelastic medium, then f IB.SP is still 
valid if appropriate Beltrami-like equations are used: see 25| for the fluid- 
elastic case. 



Step 3: high-order boundary values. The high-order boundary con- 
ditions in ( 1B.2|) and the high-order Beltrami equations (IB.SP give the under- 
determined linear systems 

\JlG\Y\ = Q, £ = 0,1. (B.6) 
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We obtain 

V,^ = K,^W,^ £ = 0,1, (B.7) 

where are the matrices built from the kernel of . The solution 

is the minimum set of independent components of the trace of U and its 
spatial derivatives up to the k-th order, on the side fl£. Injecting (1B.7P into 
the high-order jump conditions (1B.2P gives 

S^W^ = S^W^, (B.8) 

where Sf = K^. The underdetermined system (IB.Sp is solved 



W^=((S^)-^SS|Rs.) ,1, (B.9) 




where (Sj^) ^ is the least-squares pseudo-inverse of S^, Rgfc is the matrix 
containing the kernel of Sj^, and A'^ is a set of Lagrange multipliers. To build 
(S^)~^ and Rgfc, a singular value decomposition of S^^ is performed. 

Step 4: construction of modified values. We assume that U-^"^^^ 
values are known in the splitting algorithm (|30|1 . and hence the restrictions 

(2m— 2) 

Ujj are also known. Our goal is to determine a modified solution U| j at 
this time step, to be injected in the discrete operator H^. For this purpose, 
we introduce some notations. Let P be the orthogonal projection of {xj, yj) 
on r (figure E]). The coefficients of 2-D Taylor expansions around P are put 
in the matrix Ilf 



k 



(B.IO) 

where a = 0, k and /3 = 0, a; I is the 3 x 3 or 5 x 5 identity matrix, 
depending on whether (xj, yj) belongs to a fluid or a viscoelastic medium. 
By deflnition, the modifled value at (x/, yj) is 

\JIj = UIj\]1 (B.ll) 

The trace in (fRTTl) still remains to be determined in terms of the interface 

T{2n 



conditions and the numerical values \jf^ '^^ at surrounding nodes. 
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Consider the disc V is centered at P with radius q (figure [3]) • At the grid 
nodes of P fl f2o, k-ih order Taylor expansion of the solution at P, and the 
conditions flB.5|) and flB.7|) . give 



U 



■(2m-2) 



Y\k Tjfc 



^f,,.GgKgWo^ 



nf^^.GgKg (1|0) 



(B.12) 



At the grid nodes of V (1 Qi, k-th order Taylor expansion of the solution at 
P, and the interface conditions ( ]B.5|) . ( 1B.7P and ( 1B.9I) . give 



TTfc TTfe 



A^ 



(B.13) 



The equations (1B.12I) and (1B.13I) are written using an adequate matrix M 

(B.14) 



M 



V 



A" 



The radius q is chosen so that fIB.Mp is overdetermined. The least-squares 
inverse of the matrix M is denoted by M~^. Since the Lagrange multipliers 
A'^ are not involved in ( IB.llI) . is restricted to M , so that 



V 



(B.15) 



The modified value follows from flR5|) . f lRTl) . f lBTT]) and f lRlSj) : 



U 



i,j 



GS M U 



f{2m-2) 



V 



(B.16) 



V 
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